library(magrittr)
library(tidyverse)
library(Seurat)
library(readxl)
library(cowplot)
library(colorblindr)
library(viridis)
library(progeny)
library(destiny)

coi <- params$cell_type_super
cell_sort <- params$cell_sort
cell_type_major <- params$cell_type_major
louvain_resolution <- params$louvain_resolution
louvain_cluster <- params$louvain_cluster

1 Cluster markers

1.1 Major T.super markers for cell assign

### load all data ---------------------------------
source("_src/global_vars.R")

# seu_obj <- read_rds(paste0("/work/shah/isabl_data_lake/analyses/16/52/1652/celltypes/", coi, "_processed.rds"))
seu_obj <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/outs_pre/", coi, "_seurat_", louvain_resolution, ".rds"))
# seu_obj <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_highqc.rds"))

myfeatures <- c("umapharmony_1", "umapharmony_2", "sample", louvain_cluster, "doublet", "nCount_RNA", "nFeature_RNA", "percent.mt", "doublet_score")

plot_data_wrapper <- function(cluster_res) {
  cluster_res <- enquo(cluster_res)
  as_tibble(FetchData(seu_obj, myfeatures)) %>% 
    left_join(meta_tbl, by = "sample") %>% 
    rename(cluster = !!cluster_res) %>% 
    mutate(cluster = as.character(cluster),
           tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite)))
}

plot_data <- plot_data_wrapper(louvain_cluster)

1.2 Subtype currated markers

helper_f <- function(x) ifelse(is.na(x), "", x)

markers_v6_super[[coi]] %>% 
  group_by(subtype) %>% 
  mutate(rank = row_number(gene)) %>% 
  spread(subtype, gene) %>% 
  mutate_all(.funs = helper_f) %>% 
  formattable::formattable()
rank CD4.T CD4.T.CXCL13 CD4.T.KLRB1 CD4.T.reg CD8.T CD8.T.CXCL13 CD8.T.ISG Cycling.T.NK dissociated NK.CD56 NK.cytotoxic
1 CCR7 CD4 CCR6 CD4 CCL4 CCL3 IFI6 ASPM BTG1 AREG FCGR3A
2 CD4 CD40LG CD4 FOXP3 CCL5 CCL4L2 IFIT1 CENPF DNAJB1 FCER1G FGFBP2
3 CD40LG CTLA4 IL4I1 IL2RA CD8A CD8A IFIT2 HIST1H4C DUSP1 GNLY KLRD1
4 IL7R CXCL13 IL7R TNFRSF4 CD8B CRTAM IFIT3 HMGB2 EGR1 KLRC1 KLRF1
5 KLF2 FKBP5 KLRB1 TRAC CRTAM CXCL13 ISG15 MKI67 FOS KRT81 PRF1
6 LTB IL6ST LST1 CST7 FABP5 MX1 STMN1 FOSB TRDC SPON2
7 TCF7 ITM2A LTB DTHD1 GZMB MX2 TOP2A HSPA1A TYROBP
8 TPT1 MAF TNFSF13B GZMA HAVCR2 RSAD2 TUBA1B HSPA1B XCL1
9 NMB GZMH IFNG TUBB HSPA6 XCL2
10 NR3C1 GZMK LAG3 TYMS JUN
11 PDCD1 GZMM MIR155HG JUNB
12 TNFRSF4 HLA-DPB1 PHLDA1 KLF2
13 TOX2 ITM2C PTMS MT1E
14 TSHZ2 KLRG1 RBPJ MT1X
15 TRGC2 TNFRSF9

1.3 Subtype cluster markers

# marker_tbl <- read_tsv(paste0("/work/shah/isabl_data_lake/analyses/16/52/1652/celltypes/", coi, "_markers.tsv")) %>% 
#   filter(resolution == louvain_resolution)
marker_tbl <- read_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/outs_pre/", coi, "_markers_", louvain_resolution, ".tsv"))
# marker_tbl <- read_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_highqc_markers_02.tsv"))

## Hypergeometric test --------------------------------------

test_set <- marker_tbl %>% 
  group_by(cluster) %>% 
  filter(!str_detect(gene, "^RPS|^RPL")) %>% 
  slice(1:50) %>% 
  mutate(k = length(cluster)) %>% 
  ungroup %>%
  select(cluster, gene, k) %>% 
  mutate(join_helper = 1) %>% 
  group_by(cluster, join_helper, k) %>%
  nest(test_set = gene)

markers_doub_tbl <- markers_v6 %>% 
  enframe("subtype", "gene") %>% 
  filter(!(subtype %in% unique(c(coi, cell_type_major)))) %>% 
  unnest(gene) %>% 
  group_by(gene) %>% 
  filter(length(gene) == 1) %>% 
  mutate(subtype = paste0("doublet.", subtype)) %>% 
  bind_rows(tibble(subtype = "Mito.high", gene = grep("^MT-", rownames(seu_obj), value = T)))

ref_set <- markers_v6_super[[coi]] %>% 
  bind_rows(markers_doub_tbl) %>% 
  group_by(subtype) %>% 
  mutate(m = length(gene),
         n = length(rownames(seu_obj))-m,
         join_helper = 1) %>% 
  group_by(subtype, m, n, join_helper) %>%
  nest(ref_set = gene)

hyper_tbl <- test_set %>% 
  left_join(ref_set, by = "join_helper") %>% 
  group_by(cluster, subtype, m, n, k) %>%
  do(q = length(intersect(unlist(.$ref_set), unlist(.$test_set)))) %>%
  mutate(pval = 1-phyper(q = q, m = m, n = n, k = k)) %>%
  ungroup %>%
  mutate(qval = p.adjust(pval, "BH"),
         sig = qval < 0.01)

# hyper_tbl %>% 
#   group_by(subtype) %>% 
#   filter(any(qval < 0.01)) %>%
#   ggplot(aes(subtype, -log10(qval), fill = sig)) +
#   geom_bar(stat = "identity") +
#   facet_wrap(~cluster) +
#   coord_flip()
  
low_rank <- str_detect(unique(hyper_tbl$subtype), "doublet|dissociated")
subtype_lvl <- c(sort(unique(hyper_tbl$subtype)[!low_rank]), sort(unique(hyper_tbl$subtype)[low_rank]))
  
cluster_label_tbl <- hyper_tbl %>% 
  mutate(subtype = ordered(subtype, levels = subtype_lvl)) %>% 
  arrange(qval, subtype) %>%
  group_by(cluster) %>% 
  slice(1) %>% 
  mutate(subtype = ifelse(sig, as.character(subtype), paste0("unknown_", cluster))) %>% 
  select(cluster, cluster_label = subtype) %>% 
  ungroup %>% 
  mutate(cluster_label = make.unique(cluster_label, sep = "_"))

seu_obj$cluster_label <- unname(deframe(cluster_label_tbl)[as.character(unlist(seu_obj[[paste0("RNA_snn_res.", louvain_resolution)]]))])
plot_data$cluster_label <- seu_obj$cluster_label

marker_sheet <- marker_tbl %>% 
  left_join(cluster_label_tbl, by = "cluster") %>% 
  group_by(cluster_label) %>% 
  filter(!str_detect(gene, "^RPS|^RPL")) %>% 
  slice(1:50) %>% 
  mutate(rank = row_number(-avg_logFC)) %>% 
  select(cluster_label, gene, rank) %>% 
  spread(cluster_label, gene) %>% 
  mutate_all(.funs = helper_f)

formattable::formattable(marker_sheet)
rank CD4.T CD4.T.CXCL13 CD4.T.KLRB1 CD4.T.reg CD8.T CD8.T.CXCL13 CD8.T.ISG Cycling.T.NK Cycling.T.NK_1 dissociated dissociated_1 dissociated_2 doublet.Fibroblast doublet.Fibroblast_1 doublet.Monocyte doublet.Plasma.cell NK.CD56 NK.cytotoxic
1 IL7R CXCL13 KLRB1 TNFRSF4 GZMK CXCL13 IFIT3 STMN1 CENPF CCL4L2 HSPA1A KLF2 DCN IGFBP5 HLA-DRA SOX4 GNLY FGFBP2
2 CCR7 NMB IL4I1 IL2RA CD8A GZMB ISG15 TUBA1B ASPM CCL4 HSPA1B CCR7 IGFBP5 TAGLN CST3 PTCRA TYROBP FCGR3A
3 KLF2 NR3C1 IL7R FOXP3 CD8B CCL4L2 MX1 MKI67 MKI67 IFNG MT1X JUNB RBP1 ADIRF CXCL8 MAL AREG SPON2
4 EEF1B2 MAF LTB CTLA4 ITM2C MIR155HG IFIT1 HIST1H4C HMGB2 FOS DNAJB1 FOS C7 DCN SPP1 MZB1 KLRC1 PRF1
5 TPT1 FKBP5 LST1 LTB GZMH TNFRSF9 RSAD2 TUBB TOP2A FOSB MT1E DUSP1 MEG3 IGFBP7 S100A9 DNTT FCER1G KLRF1
6 EEF1A1 IL6ST TNFSF13B RTKN2 CCL5 HAVCR2 IFIT2 TOP2A UBE2C TNF HSPA6 SELL CALD1 SPARCL1 S100A8 TFDP2 TRDC GNLY
7 MAL ITM2A CCR6 BATF TRGC2 RBPJ IFI6 TYMS CCNB1 JUN HSP90AA1 IL7R IGFBP4 CALD1 HLA-DQA1 CD1E XCL1 KLRD1
8 TCF7 TSHZ2 CTSH TNFRSF18 KLRG1 LAG3 MX2 CENPF UBE2S EGR1 HSPH1 AREG RARRES2 MGP CD74 STMN1 KLRD1 NKG7
9 CD40LG CTLA4 AQP3 SAT1 GZMA IFNG IFI44L HMGB2 PTTG1 CCL3L1 HSPE1 EEF1A1 NR2F2 C11orf96 LYZ ARPP21 KRT81 CX3CR1
10 SELL CD40LG CCL20 TBC1D4 CCL4 CCL3 ISG20 ASPM TPX2 NFKBID HSPD1 CD69 SELENOP MYL9 FTL AC084033.3 XCL2 GZMB
11 GPR183 PDCD1 NFKBIA TIGIT CRTAM PTMS HERC5 NUSAP1 STMN1 CD69 HSPB1 BTG2 MDK TIMP3 APOE CDK6 IGFBP2 PLAC8
12 LDHB CD4 RORA GADD45A CST7 CD8A SAMD9L PCLAF KPNA2 NR4A2 HSP90AB1 GPR183 SOX4 MDK HLA-DQB1 MAP1A CLIC3 CLIC3
13 NOSIP LIMS1 TNFRSF25 TNFRSF1B GZMM CRTAM OAS1 HMGN2 CENPE AC020916.1 DNAJA1 PIK3IP1 ADIRF TPM2 MARCKS AC011893.1 IL2RB PLEK
14 SNHG8 TNFRSF4 CEBPD PMAIP1 HLA-DPB1 PHLDA1 TNFSF10 H2AFZ CKS2 EGR2 JUN CD55 MGP NR2F2 AIF1 GLUL CEBPD TYROBP
15 PABPC1 RNF19A TNFAIP3 UGP2 DTHD1 FABP5 STAT1 PCNA TUBA1B DUSP1 CACYBP DNAJB1 EGR1 FBLN1 C1QB ADA KRT86 PTGDS
16 NOP53 CORO1B NCR3 IKZF2 PPP1R14B TIGIT EIF2AK2 HIST1H1B HMMR TNFSF9 HSPA8 FKBP5 STAR IGF1 FTH1 AL357060.1 TXK EFHD2
17 LEF1 RBPJ SLC4A10 TNFRSF9 CD3G KRT86 OAS3 DUT DLGAP5 KLF6 FKBP4 TSC22D3 SFRP4 RAMP1 BASP1 GRASP CTSW FCER1G
18 EIF3E CPM TMIGD2 ICOS THEMIS JAML SAMD9 CLSPN CCNB2 TNFAIP3 CHORDC1 RACK1 CLU SELENOP C1QA CD1B KLRB1 CST7
19 LTB ZBED2 DPP4 LINC01943 DUSP2 CCL5 XAF1 SMC4 TUBB4B DUSP2 RGS2 LDHB ADAMTS1 RARRES2 C15orf48 MIR181A1HG MATK GZMH
20 RACK1 AC004585.1 TPT1 IL32 CD3D LINC01871 GBP1 ATAD2 NUSAP1 BTG2 FOS SARAF TIMP2 LUM FN1 CCDC26 CCL3 ADGRG1
21 NACA TOX2 MYBL1 SOX4 LYAR CXCR6 EPSTI1 SMC2 BIRC5 IER2 ANXA1 CXCR4 TCEAL4 IGFBP6 MNDA JCHAIN CD7 CCL3
22 UBA52 DUSP4 S100A4 ARID5B TC2N TNIP3 IFI44 TPX2 HMGN2 PPP1R15A DNAJB4 PLAC8 C1R CARMN G0S2 VIPR2 NKG7 HOPX
23 TOMM7 AHI1 CD40LG CD27 EOMES PDCD1 PLSCR1 TMPO ARL6IP1 NR4A1 PPP1R15A EEF1B2 WFDC2 IGFBP4 APOC1 ID1 CD63 IGFBP7
24 SOCS3 ICA1 SPOCK2 BIRC3 CXCR6 HLA-DRB1 IFI35 HELLS CDC20 NFKBIA ZFAND2A TCF7 FHL2 DST SOD2 SOCS2 HOPX ZEB2
25 JUNB ARID5B LINC01871 LAYN HLA-DPA1 GAPDH USP18 UBE2C CDKN3 JUND FOSB TPT1 IFITM3 CAV1 NPC2 RCAN1 TMIGD2 PRSS23
26 SERINC5 CD84 ERN1 CORO1B HLA-DRB1 FAM3C OASL NASP SMC4 GADD45B DNAJA4 FOSB PEG3 ID4 LST1 CLDN5 CMC1 AKR1C3
27 TMEM123 CCDC50 JAML TYMP SLF1 CTLA4 MT2A DEK CKS1B TSC22D3 DUSP1 NOSIP C1S NUPR1 GSN CASC15 TNFRSF18 CD247
28 EEF2 IGFL2 FKBP11 DUSP4 APOBEC3G SPRY1 HELZ2 RRM2 KIF20B TAGAP TSPYL2 SC5D LUM CSRP2 MS4A6A PFKFB2 KLRC2 MYBL1
29 FXYD5 RGS1 IFNGR1 CD4 PPP2R5C CCND2 TRIM22 CKS1B GTSE1 ZFP36 SERPINH1 NACA SERPINF1 CDKN1C IL1B GALNT2 SRGAP3 AREG
30 TSHZ2 BATF MGAT4A ENTPD1 KIAA1551 CD63 CMPK2 KNL1 TUBA1C ZFP36L1 UBC AP3M2 CEBPD PGR PSAP HES4 GSTP1 C1orf21
31 TRABD2A SRGN S100A6 CTSC F2R CD8B LY6E HMGB1 TUBB RGCC UBB BTG1 AKAP12 COL6A1 GRN MARCKSL1 LAT2 S1PR5
32 ANK3 CH25H ELK3 MIR4435-2HG SLAMF7 GZMH PARP14 MCM7 SGO2 IL7R AHSA1 NOP53 TIMP1 COL6A2 CD83 TP53INP1 GZMB CTSW
33 SARAF SPOCK2 EEF1A1 LINC02099 CXCR4 ENTPD1 NT5C3A FABP5 H2AFZ CRTAM NEU1 ZBTB16 CST3 EMX2 CTSH APBA2 LINC00996 ABHD17A
34 AQP3 ZNRF1 KIT MAGEH1 SH2D1A SRGAP3 DDX58 UBE2S JPT1 KDM6B DNAJB6 ERAP2 NR2F1 PALLD GLUL TSHR PRF1 KLF2
35 RIPOR2 CHN1 LTC4S SPOCK2 FAM102A VCAM1 IRF7 GAPDH CEP55 NR4A3 GADD45B CCND3 DLK1 PPP1R14A MEF2C NREP KLRF1 PTPN12
36 AP3M2 TNFRSF25 RUNX2 PHACTR2 CD52 ID2 RNF213 CXCL13 NUF2 ATF3 KLF6 EEF1D SERPING1 RBP1 CYBB MME NCAM1 TTC38
37 ZFAS1 CD200 RORC CARD16 YBX3 HLA-DRA DDX60L RANBP1 KIF14 DDX3X BTG2 PPP1R15A BEX3 C7 CTSB AC002454.1 CXXC5 KLRB1
38 LINC02273 METTL8 ZBTB16 S100A4 STK17A ITGAE DDX60 H2AFX KNL1 DUSP6 JUNB PLK3 C11orf96 MFGE8 CSF3R CHI3L2 SH2D1B CCL4
39 EIF4B RILPL2 FAM241A STAM CCR5 DUSP4 SAT1 MCM3 CENPA GZMK TNF ZFP36L2 FILIP1L KANK2 CD14 SMIM3 MCTP2 PTGDR
40 TOB1 TNFRSF18 PDE4D GLRX GPR174 LYST PPM1K EZH2 CDK1 KLRG1 ERN1 EEF2 RARRES1 NR2F1 C1QC SSBP2 IFITM3 ITGB2
41 SESN3 SLA IL23R SPATS2L COTL1 TNFSF4 PNPT1 TUBB4B HMGB3 JUNB JUND HNRNPA1 COL1A2 MFAP4 SGK1 UHRF1 ZNF683 XBP1
42 NSA2 SMCO4 B3GALT2 AC005224.3 CCL4L2 NDFIP2 PARP9 H2AFV CDCA8 RASGEF1B ATF3 VSIR NUPR1 COL1A1 SPI1 LRRC28 ITGA1 CEP78
43 PASK BTLA CERK MAF CD3E AKAP5 IFIH1 CDK1 PLK1 ANXA1 DEDD2 PASK TCEAL9 PBX1 FCGRT BCL11A IFITM2 ARL4C
44 TNFRSF25 NAP1L4 TLE1 PBXIP1 TUBA4A GOLIM4 SP110 HIST1H1D AURKA MCL1 CD69 EIF3H MARCKSL1 GSN EGR1 SCAI CCL5 BIN2
45 FAU FYB1 EEF1B2 F5 PECAM1 CD27 OAS2 HNRNPAB TROAP CXCR4 CLK1 TXNIP SLC40A1 PDGFRB FCGR2A ATP6AP1L ITGAX LITAF
46 EEF1D MIR155HG CFH SLAMF1 ARAP2 SNAP47 C19orf66 CENPE H2AFV IER5 IER5L LDLRAP1 CFH SERPINF1 PLAUR RUFY3 CD38 TRDC
47 LDLRAP1 PTPN13 PERP BTG3 ITGA1 RGS1 STAT2 TK1 KIF2C MYADM H3F3B CMTM8 APOE SFRP1 CPVL CD79A SAMD3 TXK
48 CTSL SESN3 PLAT TRAC JAML HLA-DPA1 LAG3 ZWINT MAD2L1 CD8A NR4A1 SCML1 GNG11 LHFPL6 MS4A7 GNA15 SLC16A3 MYOM2
49 ITGA6 BIRC3 KIF5C IL1R1 CD84 LINC02446 LAP3 BIRC5 NUCKS1 PTGER4 CXCR4 LINC00402 CDKN1C PLAC9 ALDH2 HHIP-AS1 CAPN12 GZMM
50 PFDN5 TP53INP1 PLCB1 DNPH1 AOAH SAMSN1 APOL6 PTTG1 RAD21 ZFP36L2 EIF4A2 RIPK2 NBL1 SERPING1 SERPINA1 GSTM3 CD247 CD300A
write_tsv(marker_sheet, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_sheet.tsv"))

2 Clusters

2.1 sizes

enframe(sort(table(seu_obj$cluster_label))) %>% 
  mutate(name = ordered(name, levels = rev(name))) %>% 
  ggplot() +
  geom_bar(aes(name, value), stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(y = c("#cells"), x = "cluster")

2.2 UMAP

alpha_lvl <- ifelse(nrow(plot_data) < 20000, 0.2, 0.1)
pt_size <- ifelse(nrow(plot_data) < 20000, 0.2, 0.05)

common_layers_disc <- list(  
  geom_point(size = pt_size, alpha = alpha_lvl),
  NoAxes(),
  guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))),
  labs(color = "")
)

common_layers_cont <- list(  
  geom_point(size = pt_size, alpha = alpha_lvl),
  NoAxes(),
  scale_color_gradientn(colors = viridis(9)),
  guides(color = guide_colorbar())
)

ggplot(plot_data, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
  common_layers_disc +
  #facet_wrap(~therapy) +
  ggtitle("Sub cluster")

3 Filtering out doublet clusters

3.1 UMAP

my_subtypes <- names(clrs$cluster_label[[coi]])
my_subtypes <- c(my_subtypes, unlist(lapply(paste0("_", 1:3), function(x) paste0(my_subtypes, x)))) %>% .[!str_detect(., "doublet|dissociated")]
my_subtypes <- my_subtypes[my_subtypes %in% unique(seu_obj$cluster_label)]
my_subtypes <- my_subtypes[my_subtypes %in% names(clrs$cluster_label[[coi]])]

cells_to_keep <- colnames(seu_obj)[seu_obj$cluster_label %in% my_subtypes]
seu_obj_sub <- subset(seu_obj, cells = cells_to_keep)
seu_obj_sub <- RunUMAP(seu_obj_sub, dims = 1:50, reduction = "harmony", reduction.name = "umapharmony", reduction.key = "umapharmony_")
seu_obj_sub$cluster_label <- seu_obj$cluster_label[colnames(seu_obj) %in% colnames(seu_obj_sub)]
write_rds(seu_obj_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered.rds"))
seu_obj_sub <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered.rds"))

plot_data_sub <- as_tibble(FetchData(seu_obj_sub, c(myfeatures, "cluster_label"))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
  mutate(cell_id = colnames(seu_obj_sub),
         cluster_label = ordered(cluster_label, levels = my_subtypes),
         )
  
if (cell_sort == "CD45+") {
plot_data_sub <- filter(plot_data_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
}

if (cell_sort == "CD45-") {
plot_data_sub <- filter(plot_data_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
}

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
  common_layers_disc +
  ggtitle("Sub cluster") +
  #facet_wrap(~cluster_label) +
  scale_color_manual(values = clrs$cluster_label[[coi]])

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = patient_id_short)) + 
  common_layers_disc +
  ggtitle("Patient") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$patient_id_short)

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
  # geom_point(aes(umapharmony_1, umapharmony_2), 
  #            color = "grey90", size = 0.01, 
  #            data = select(plot_data_sub, -tumor_supersite)) +
  common_layers_disc +
  ggtitle("Site") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$tumor_supersite)

write_tsv(select(plot_data_sub, cell_id, everything(), -umapharmony_1, -umapharmony_2, -contains("RNA_")), paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_embedding.tsv"))

4 subcluster big clusters

# big_clusters <- list(
#   NK = c("NK"),
#   CD8.T = c("CD8.T"),
#   CD4.T = c("CD4.T")
# )
# 
# cells_to_keep_list <- lapply(big_clusters, function(x) colnames(seu_obj_sub)[seu_obj_sub$cluster_label %in% x])
# 
# # seu_list <- lapply(cells_to_keep_list, function(x) subset(seu_obj_sub, cells = x))
# # 
# # preprocess_wrapper <- . %>%
# #   NormalizeData() %>%
# #   FindVariableFeatures() %>%
# #   ScaleData() %>%
# #   RunPCA(verbose = F) %>%
# #   RunHarmony(group.by.vars = "patient_id") %>%
# #   FindNeighbors(reduction = "harmony", dims = 1:50) %>%
# #   FindClusters(res = 0.1)
# # 
# # seu_list <- lapply(seu_list, preprocess_wrapper)
# # write_rds(seu_list, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub_clusters.rds"))
# 
# seu_list <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub_clusters.rds"))
# 
# myfeatures <- c("umapharmony_1", "umapharmony_2", "sample", "RNA_snn_res.0.1")
# 
# plot_data_clusters <- lapply(seu_list, function(x) as_tibble(FetchData(x, myfeatures))) %>% 
#   bind_rows(.id = "big_cluster")
# 
# ggplot(plot_data_clusters, aes(umapharmony_1, umapharmony_2, 
#                                color = RNA_snn_res.0.1)) + 
#   common_layers_disc +
#   facet_grid(~big_cluster)

4.1 add module scores and pathway scores

signature_modules <- read_excel("_data/small/signatures/SPECTRUM v5 sub cluster markers.xlsx", sheet = 2, skip = 1, range = "M2:P100") %>% 
  gather(module, gene) %>% 
  na.omit() %>% 
  group_by(module) %>% 
  do(gene = c(.$gene)) %>% 
  {setNames(.$gene, .$module)}

signature_modules$ISG.module <- c("CCL5", "CXCL10", "IFNA1", "IFNB1", "ISG15", "IFI27L2", "SAMD9L")

## compute expression module scores
for (i in 1:length(signature_modules)) {
  seu_obj_sub <- AddModuleScore(seu_obj_sub, features = signature_modules[i], name = names(signature_modules)[i])
  seu_obj_sub[[names(signature_modules)[i]]] <- seu_obj_sub[[paste0(names(signature_modules)[i], "1")]]
  seu_obj_sub[[paste0(names(signature_modules)[i], "1")]] <- NULL
  print(paste(names(signature_modules)[i], "DONE"))
}
## [1] "CD8.Cytotoxic DONE"
## [1] "CD8.Dysfunctional DONE"
## [1] "CD8.Naive DONE"
## [1] "CD8.Predysfunctional DONE"
## [1] "ISG.module DONE"
## compute progeny scores
progeny_list <- seu_obj_sub@assays$RNA@data[VariableFeatures(seu_obj_sub),] %>% 
  as.matrix %>% 
  progeny %>% 
  as.data.frame %>% 
  as.list

names(progeny_list) <- make.names(paste0(names(progeny_list), ".pathway"))

for (i in 1:length(progeny_list)) {
  seu_obj_sub <- AddMetaData(seu_obj_sub, metadata = progeny_list[[i]], 
                             col.name = names(progeny_list)[i])
}

write_rds(seu_obj_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_annotated.rds"))

4.2 marker heatmap

marker_top_tbl <- marker_sheet[,-1] %>% 
  slice(1:10) %>% 
  as.list %>% 
  .[!str_detect(names(.), "doublet")] %>% 
  enframe("cluster_label_x", "gene") %>% 
  unnest(gene)

plot_data_markers <- as_tibble(FetchData(seu_obj_sub, c("cluster_label", myfeatures, unique(marker_top_tbl$gene)))) %>% 
  gather(gene, value, -c(1:(length(myfeatures)+1))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
  mutate(cluster_label = ordered(cluster_label, levels = my_subtypes)) %>% 
  group_by(cluster_label, gene) %>% 
  summarise(value = mean(value, na.rm = T)) %>% 
  group_by(gene) %>% 
  mutate(value = scales::rescale(value)) %>% 
  left_join(marker_top_tbl, by = "gene") %>% 
  mutate(cluster_label_x = ordered(cluster_label_x, levels = rev(names(clrs$cluster_label[[coi]]))))

ggplot(plot_data_markers) +
  geom_tile(aes(gene, cluster_label, fill = value)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  facet_grid(~cluster_label_x, scales = "free", space = "free") +
  scale_fill_gradientn(colors = viridis(9)) +
  labs(fill = "Scaled\nexpression") +
  theme(aspect.ratio = 1,
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank())

# ggsave(paste0("_fig/002_marker_heatmap_", coi, ".pdf"), width = nrow(marker_top_tbl)/6, height = 5)

4.3 composition

4.3.1 per site

comp_site_tbl <- plot_data_sub %>%
  filter(!is.na(tumor_supersite)) %>% 
  group_by(cluster_label, tumor_supersite) %>%
  tally %>%
  group_by(tumor_supersite) %>%
  mutate(nrel = n/sum(n)*100) %>%
  ungroup

pnrel_site <- ggplot(comp_site_tbl) +
  geom_bar(aes(tumor_supersite, nrel, fill = cluster_label),
           stat = "identity", position = position_stack()) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(fill = "Cluster", y = "Fraction [%]", x = "") +
  scale_fill_manual(values = clrs$cluster_label[[coi]])

pnabs_site <- ggplot(comp_site_tbl) +
  geom_bar(aes(tumor_supersite, n, fill = cluster_label),
           stat = "identity", position = position_stack()) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(fill = "Cluster", y = "# cells", x = "") +
  scale_fill_manual(values = clrs$cluster_label[[coi]])

plot_grid(pnabs_site, pnrel_site, ncol = 2, align = "h")

# ggsave(paste0("_fig/02_deep_dive_", coi, "_comp_site.pdf"), width = 8, height = 4)

4.3.2 per sample

comp_tbl_sample_sort <- plot_data_sub %>% 
  group_by(tumor_subsite, tumor_supersite, patient_id, consensus_signature, therapy, cluster_label) %>% 
  tally %>% 
  group_by(tumor_subsite, tumor_supersite, patient_id, consensus_signature, therapy) %>% 
  mutate(nrel = n/sum(n)*100,
         nsum = sum(n),
         log10n = log10(n),
         label_supersite = "Site",
         label_mutsig = "Signature",
         label_therapy = "Rx") %>% 
  ungroup %>% 
  arrange(desc(therapy), tumor_supersite) %>% 
  mutate(tumor_subsite_rx = paste0(tumor_subsite, "_", therapy)) %>% 
  mutate(tumor_subsite = ordered(tumor_subsite, levels = unique(tumor_subsite)),
         tumor_subsite_rx = ordered(tumor_subsite_rx, levels = unique(tumor_subsite_rx))) %>% 
  arrange(patient_id) %>% 
  mutate(label_patient_id = ifelse(as.logical(as.numeric(fct_inorder(as.character(patient_id)))%%2), "Patient1", "Patient2"))

sample_id_x_tbl <- plot_data_sub %>% 
  mutate(sort_short_x = cell_sort) %>% 
  distinct(patient_id, sort_short_x, tumor_subsite, therapy, sample) %>% 
  unite(sample_id_x, patient_id, sort_short_x, tumor_subsite, therapy) %>% 
  arrange(sample_id_x)

comp_tbl_sample_sort %>% 
  mutate(sort_short_x = cell_sort) %>% 
  unite(sample_id_x, patient_id, sort_short_x, tumor_subsite_rx) %>% 
  select(sample_id_x, cluster_label, n, nrel, nsum) %>% 
  left_join(sample_id_x_tbl, by = "sample_id_x") %>% 
  write_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_subtype_compositions.tsv"))
ybreaks <- c(500, 1000, 2000, 4000, 6000, 8000, 10000, 15000, 20000)

max_cells_per_sample <- max(comp_tbl_sample_sort$nsum)
ymaxn <- ybreaks[max_cells_per_sample < ybreaks][1]

comp_plot_wrapper <- function(y = "nrel", switch = NULL) {
  if (y == "nrel") ylab <- paste0("Fraction\nof cells [%]")
  if (y == "n") ylab <- paste0("Number\nof cells")
  p <- ggplot(comp_tbl_sample_sort, 
              aes_string("tumor_subsite_rx", y, fill = "cluster_label")) + 
    facet_grid(~patient_id, space = "free", scales = "free", switch = switch) +
    coord_cartesian(clip = "off") + 
    scale_fill_manual(values = clrs$cluster_label[[coi]]) + 
    theme(axis.text.x = element_blank(),
          axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5, 
                                      margin = margin(0, -0.4, 0, 0, unit = "npc")),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          axis.line.x = element_blank(),
          strip.text.y = element_blank(),
          strip.text.x = element_blank(),
          strip.background.y = element_blank(),
          strip.background.x = element_blank(),
          plot.margin = margin(0, 0, 0, 0, "npc")) + 
    labs(x = "",
         y = ylab) +
    guides(fill = FALSE)
  if (y == "nrel") p <- p + 
    geom_bar(stat = "identity") +
    scale_y_continuous(expand = c(0, 0), 
                       breaks = c(0, 50, 100), 
                       labels = c("0", "50", "100"))
  if (y == "n") p <- p + 
    geom_bar(stat = "identity", position = position_stack()) +
    scale_y_continuous(expand = c(0, 0), 
                       breaks = c(0, ymaxn/2, ymaxn),
                       limits = c(0, ymaxn),
                       labels = c("", ymaxn/2, ymaxn)) +
    expand_limits(y = c(0, ymaxn)) +
    theme(panel.grid.major.y = element_line(linetype = 1, color = "grey90", size = 0.5))
  return(p)
} 

common_label_layers <- list(
  geom_tile(color = "white", size = 0.15),
  theme(axis.text.x = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.line.x = element_blank(),
        strip.text = element_blank(),
        strip.background = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "npc")),
  scale_y_discrete(expand = c(0, 0)),
  labs(y = ""),
  guides(fill = FALSE),
  facet_grid(~patient_id, 
             space = "free", scales = "free")
)

comp_label_site <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_supersite, patient_id), 
       aes(tumor_subsite_rx, label_supersite, 
           fill = tumor_supersite)) + 
  scale_fill_manual(values = clrs$tumor_supersite) +
  common_label_layers

comp_label_rx <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_therapy, consensus_signature, patient_id), 
       aes(tumor_subsite_rx, label_therapy, 
           fill = therapy)) + 
  scale_fill_manual(values = c(`post-Rx` = "gold3", `pre-Rx` = "steelblue")) +
  common_label_layers

comp_label_mutsig <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_mutsig, consensus_signature, patient_id), 
       aes(tumor_subsite_rx, label_mutsig, 
           fill = consensus_signature)) + 
  scale_fill_manual(values = clrs$consensus_signature) +
  common_label_layers

patient_label_tbl <- distinct(comp_tbl_sample_sort, patient_id, .keep_all = T)

comp_label_patient_id <- ggplot(patient_label_tbl, aes(tumor_subsite_rx, label_patient_id)) + 
  scale_fill_manual(values = clrs$consensus_signature) +
  geom_text(aes(tumor_subsite_rx, label_patient_id, label = patient_id)) +
  facet_grid(~patient_id, 
             space = "free", scales = "free") +
  coord_cartesian(clip = "off") + 
  theme_void() +
  theme(strip.text = element_blank(),
        strip.background = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "npc"))

hist_plot_wrapper <- function(x = "nrel") {
  if (x == "nrel") {
    xlab <- paste0("Fraction of cells [%]")
    bw <- 5
  }
  if (x == "log10n") {
    xlab <- paste0("Number of cells")
    bw <- 0.2
  }
  p <- ggplot(comp_tbl_sample_sort) +
    ggridges::geom_density_ridges(
      aes_string(x, "cluster_label", fill = "cluster_label"), color = "black",
      stat = "binline", binwidth = bw, scale = 3) +
    facet_grid(label_supersite~., 
               space = "free", scales = "free") +
    scale_fill_manual(values = clrs$cluster_label[[coi]]) + 
    theme(axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.title.y = element_blank(),
          axis.line.y = element_blank(),
          strip.text = element_blank(),
          strip.background = element_blank(),
          plot.margin = margin(0, 0, 0, 0, "npc")) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_discrete(expand = c(0, 0)) +
    guides(fill = FALSE) +
    labs(x = xlab)
  if (x == "log10n") p <- p + expand_limits(x = c(0, 3)) + 
    scale_x_continuous(expand = c(0, 0), 
                       labels = function(x) parse(text = paste("10^", x)))
  return(p)
}

pcomp1 <- comp_plot_wrapper("n")
pcomp2 <- comp_plot_wrapper("nrel")

pcomp_grid <- plot_grid(comp_label_patient_id, 
                        pcomp1, pcomp2, 
                        comp_label_site, comp_label_rx, comp_label_mutsig,
                        ncol = 1, align = "v", 
                        rel_heights = c(0.15, 0.33, 0.33, 0.06, 0.06, 0.06))

phist1 <- hist_plot_wrapper("log10n")

pcomp_hist_grid <- ggdraw() +
  draw_plot(pcomp_grid, x = 0.01, y = 0, width = 0.85, height = 1) +
  draw_plot(phist1, x = 0.87, y = 0.05, width = 0.12, height = 0.8)

pcomp_hist_grid

# ggsave(paste0("_fig/02_composition_v6_",coi,".pdf"), pcomp_hist_grid, width = 10, height = 2)

4.3.3 site specific cluster enrichment

comp_tbl_z <- comp_tbl_sample_sort %>% 
  filter(therapy == "pre-Rx",
         !(tumor_supersite %in% c("Ascites", "Other"))) %>% 
  group_by(patient_id, cluster_label) %>% 
  arrange(patient_id, cluster_label, nrel) %>% 
  mutate(rank = row_number(nrel),
         z_rank = scales::rescale(rank)) %>% 
  mutate(mean_nrel = mean(nrel, na.rm = T),
         sd_nrel = sd(nrel, na.rm = T),
         z_nrel = (nrel - mean_nrel) / sd_nrel) %>% 
  ungroup()

ggplot(comp_tbl_z) +
  geom_boxplot(aes(tumor_supersite, z_nrel, color = tumor_supersite), 
               outlier.shape = NA) +
  geom_point(aes(tumor_supersite, z_nrel, color = tumor_supersite), position = position_jitter(), size = 0.1) +
  facet_grid(~cluster_label, scales = "free_x", space = "free_x") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        aspect.ratio = 5) +
  scale_color_manual(values = clrs$tumor_supersite)

ggplot(comp_tbl_z) +
  geom_boxplot(aes(tumor_supersite, z_rank, color = tumor_supersite), 
               outlier.shape = NA) +
  geom_point(aes(tumor_supersite, z_rank, color = tumor_supersite), position = position_jitter(), size = 0.1) +
  facet_grid(~cluster_label, scales = "free_x", space = "free_x") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        aspect.ratio = 5) +
  scale_color_manual(values = clrs$tumor_supersite)

4.4 sub cluster CD8 cells

# seu_obj_cd8 <- seu_obj_sub %>%
#   subset(subset = cluster_label %in% c("CD8.T")) %>%
#   FindNeighbors(reduction = "harmony", dims = 1:50) %>%
#   FindClusters(resolution = 0.2) %>%
#   RunUMAP(dims = 1:50, reduction = "harmony", reduction.name = "umapharmony", reduction.key = "umapharmony_")

seu_obj_cd8 <- seu_obj_sub %>%
  subset(subset = cluster_label %in% my_subtypes[str_detect(my_subtypes, "CD8.T")]) %>% 
  RunUMAP(dims = 1:50, reduction = "harmony", reduction.name = "umapharmony", reduction.key = "umapharmony_")

write_rds(seu_obj_cd8, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed.rds")
seu_obj_cd8 <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed.rds")

# marker_tbl <- FindAllMarkers(seu_obj_cd8, only.pos = T)
# write_tsv(as_tibble(marker_tbl), "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_markers.tsv")
# marker_tbl <- read_tsv("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_markers.tsv")
# 
# ## Hypergeometric test --------------------------------------
# 
# test_set <- marker_tbl %>% 
#   group_by(cluster) %>% 
#   filter(!str_detect(gene, "^RPS|^RPL")) %>% 
#   slice(1:50) %>% 
#   mutate(k = length(cluster)) %>% 
#   ungroup %>%
#   select(cluster, gene, k) %>% 
#   mutate(join_helper = 1) %>% 
#   group_by(cluster, join_helper, k) %>%
#   nest(test_set = gene)
# 
# markers_doub_tbl <- markers_v6 %>% 
#   enframe("subtype", "gene") %>% 
#   filter(!(subtype %in% unique(c(coi, cell_type_major)))) %>% 
#   unnest(gene) %>% 
#   group_by(gene) %>% 
#   filter(length(gene) == 1) %>% 
#   mutate(subtype = paste0("doublet.", subtype)) %>% 
#   bind_rows(tibble(subtype = "Mito.high", gene = grep("^MT-", rownames(seu_obj), value = T)))
# 
# ref_set <- markers_v6_super[["CD8.T"]] %>% 
#   bind_rows(markers_doub_tbl) %>% 
#   group_by(subtype) %>% 
#   mutate(m = length(gene),
#          n = length(rownames(seu_obj))-m,
#          join_helper = 1) %>% 
#   group_by(subtype, m, n, join_helper) %>%
#   nest(ref_set = gene)
# 
# hyper_tbl <- test_set %>% 
#   left_join(ref_set, by = "join_helper") %>% 
#   group_by(cluster, subtype, m, n, k) %>%
#   do(q = length(intersect(unlist(.$ref_set), unlist(.$test_set)))) %>%
#   mutate(pval = 1-phyper(q = q, m = m, n = n, k = k)) %>%
#   ungroup %>%
#   mutate(qval = p.adjust(pval, "BH"),
#          sig = qval < 0.01)
# 
# # hyper_tbl %>% 
# #   group_by(subtype) %>% 
# #   filter(any(qval < 0.01)) %>%
# #   ggplot(aes(subtype, -log10(qval), fill = sig)) +
# #   geom_bar(stat = "identity") +
# #   facet_wrap(~cluster) +
# #   coord_flip()
#   
# low_rank <- str_detect(unique(hyper_tbl$subtype), "Mito|doublet")
# subtype_lvl <- c(sort(unique(hyper_tbl$subtype)[!low_rank]), sort(unique(hyper_tbl$subtype)[low_rank]))
#   
# cluster_label_tbl <- hyper_tbl %>% 
#   mutate(subtype = ordered(subtype, levels = subtype_lvl)) %>% 
#   arrange(qval, subtype) %>%
#   group_by(cluster) %>% 
#   slice(1) %>% 
#   mutate(subtype = ifelse(sig, as.character(subtype), paste0("unknown_", cluster))) %>% 
#   select(cluster, cluster_label = subtype) %>% 
#   ungroup %>% 
#   mutate(cluster_label = make.unique(cluster_label, sep = "_"))
# 
# marker_sheet <- marker_tbl %>% 
#   left_join(cluster_label_tbl, by = "cluster") %>% 
#   group_by(cluster_label) %>% 
#   filter(!str_detect(gene, "^RPS|^RPL")) %>% 
#   slice(1:50) %>% 
#   mutate(rank = row_number()) %>% 
#   select(cluster_label, gene, rank) %>% 
#   spread(cluster_label, gene) %>% 
#   mutate_all(.funs = helper_f)
# 
# formattable::formattable(marker_sheet)
# write_tsv(marker_sheet, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T.cell_marker_sheet.tsv"))
# seu_obj_cd8$cluster_label <- unname(deframe(cluster_label_tbl)[as.character(unlist(seu_obj_cd8[[paste0("RNA_snn_res.", louvain_resolution)]]))])
# 
# seu_obj_sub_sub <- seu_obj_sub
# 
# cluster_label_tbl_1 <- as_tibble(cbind(cell_id=colnames(seu_obj_sub), FetchData(seu_obj_sub, c("cluster_label"))))
# cluster_label_tbl_2 <- as_tibble(cbind(cell_id=colnames(seu_obj_cd8), FetchData(seu_obj_cd8, c("cluster_label"))))
# 
# cluster_label_tbl <- left_join(cluster_label_tbl_1, cluster_label_tbl_2, by = "cell_id") %>%
#   mutate(cluster_label = ifelse(is.na(cluster_label.y), cluster_label.x, cluster_label.y))
# 
# seu_obj_sub_sub$cluster_label <- cluster_label_tbl$cluster_label
# 
# seu_obj_sub_sub <- subset(seu_obj_sub_sub, subset = cluster_label != "doublet.Plasma.cell")
# write_rds(seu_obj_sub_sub, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/T.cell_processed_filtered_sub.rds")
root_cell <- "SPECTRUM-OV-036_S1_CD45P_PELVIC_PERITONEUM_CGGACACCAACGACTT"
seu_obj_cd8_sub <- subset(seu_obj_cd8, subset = cluster_label != "doublet.Plasma.cell" & cluster_label != "gd.T.cell")
seu_obj_cd8_sub <- subset(seu_obj_cd8_sub, cells = c(root_cell, colnames(seu_obj_cd8_sub)[colnames(seu_obj_cd8_sub)!=root_cell][-1]))

dc_obj <- DiffusionMap(seu_obj_cd8_sub@reductions$harmony@cell.embeddings, k = 100)
dc_mat <- dc_obj@eigenvectors
colnames(dc_mat) <- paste0("DC_", 1:ncol(dc_mat))
seu_obj_cd8_sub[["DC"]] <- CreateDimReducObject(embeddings = dc_mat, key = "DC_", assay = DefaultAssay(seu_obj_cd8_sub))

write_rds(seu_obj_cd8_sub, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed_filtered.rds")

dpt_obj <- DPT(dc_obj)
write_rds(dpt_obj, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_sub_dpt.rds")
seu_obj_cd8_sub$DPT1 <- dpt_obj$DPT1

write_rds(seu_obj_cd8_sub, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed_filtered.rds")
seu_obj_cd8_sub <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed_filtered.rds")
# plot_data_sub_sub <- as_tibble(FetchData(seu_obj_sub_sub, c(myfeatures, "cluster_label"))) %>% 
#   left_join(meta_tbl, by = "sample") %>% 
#   mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
#          tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
#   mutate(cell_id = colnames(seu_obj_sub_sub),
#          cluster_label = ordered(cluster_label, levels = my_subtypes),
#          )
#   
# if (cell_sort == "CD45+") {
# plot_data_sub_sub <- filter(plot_data_sub_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
# }
# 
# if (cell_sort == "CD45-") {
# plot_data_sub_sub <- filter(plot_data_sub_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
# }
# 
# ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
#   common_layers_disc +
#   ggtitle("Sub cluster") +
#   #facet_wrap(~cluster_label) +
#   scale_color_manual(values = clrs$cluster_label[[coi]])
# 
# ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = patient_id_short)) + 
#   common_layers_disc +
#   ggtitle("Patient") +
#   # facet_wrap(~therapy) +
#   scale_color_manual(values = clrs$patient_id_short)
# 
# ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
#   # geom_point(aes(umapharmony_1, umapharmony_2), 
#   #            color = "grey90", size = 0.01, 
#   #            data = select(plot_data_sub_sub, -tumor_supersite)) +
#   common_layers_disc +
#   ggtitle("Site") +
#   # facet_wrap(~therapy) +
#   scale_color_manual(values = clrs$tumor_supersite)
# 
# write_tsv(select(plot_data_sub_sub, cell_id, everything(), -umapharmony_1, -umapharmony_2, -contains("RNA_")), paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_embedding_sub.tsv"))

4.4.1 UMAP & DC

plot_data_cd8_sub <- as_tibble(FetchData(seu_obj_cd8_sub, c(myfeatures, "cluster_label", "DC_1", "DC_2"))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
  mutate(cell_id = colnames(seu_obj_cd8_sub),
         cluster_label = ordered(cluster_label, levels = my_subtypes),
         )
  
if (cell_sort == "CD45+") {
plot_data_cd8_sub <- filter(plot_data_cd8_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
}

if (cell_sort == "CD45-") {
plot_data_cd8_sub <- filter(plot_data_cd8_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
}

ggplot(plot_data_cd8_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
  common_layers_disc +
  ggtitle("Sub cluster") +
  #facet_wrap(~cluster_label) +
  scale_color_manual(values = clrs$cluster_label[[coi]])

ggplot(plot_data_cd8_sub, aes(umapharmony_1, umapharmony_2, color = patient_id_short)) + 
  common_layers_disc +
  ggtitle("Patient") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$patient_id_short)

ggplot(plot_data_cd8_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
  common_layers_disc +
  ggtitle("Site") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$tumor_supersite)

ggplot(plot_data_cd8_sub, aes(DC_1, DC_2, color = cluster_label)) + 
  common_layers_disc +
  ggtitle("Sub cluster") +
  #facet_wrap(~cluster_label) +
  scale_color_manual(values = clrs$cluster_label[[coi]])

ggplot(plot_data_cd8_sub, aes(DC_1, DC_2, color = patient_id_short)) + 
  common_layers_disc +
  ggtitle("Patient") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$patient_id_short)

ggplot(plot_data_cd8_sub, aes(DC_1, DC_2, color = tumor_supersite)) + 
  common_layers_disc +
  ggtitle("Site") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$tumor_supersite)

5 session info

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.2 (2019-12-12)
##  os       Debian GNU/Linux 10 (buster)
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       Etc/UTC                     
##  date     2021-01-20                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package              * version    date       lib
##  abind                  1.4-5      2016-07-21 [2]
##  ape                    5.3        2019-03-17 [2]
##  assertthat             0.2.1      2019-03-21 [2]
##  backports              1.1.10     2020-09-15 [1]
##  bibtex                 0.4.2.2    2020-01-02 [2]
##  Biobase                2.46.0     2019-10-29 [2]
##  BiocGenerics           0.32.0     2019-10-29 [2]
##  BiocParallel           1.20.1     2019-12-21 [2]
##  bitops                 1.0-6      2013-08-17 [2]
##  boot                   1.3-24     2019-12-20 [3]
##  broom                  0.7.2      2020-10-20 [1]
##  callr                  3.4.2      2020-02-12 [1]
##  car                    3.0-8      2020-05-21 [1]
##  carData                3.0-4      2020-05-22 [1]
##  caTools                1.17.1.4   2020-01-13 [2]
##  cellranger             1.1.0      2016-07-27 [2]
##  class                  7.3-15     2019-01-01 [3]
##  cli                    2.0.2      2020-02-28 [1]
##  cluster                2.1.0      2019-06-19 [3]
##  codetools              0.2-16     2018-12-24 [3]
##  colorblindr          * 0.1.0      2020-01-13 [2]
##  colorspace           * 1.4-2      2019-12-29 [2]
##  cowplot              * 1.0.0      2019-07-11 [2]
##  crayon                 1.3.4      2017-09-16 [1]
##  curl                   4.3        2019-12-02 [2]
##  data.table             1.12.8     2019-12-09 [2]
##  DBI                    1.1.0      2019-12-15 [2]
##  dbplyr                 2.0.0      2020-11-03 [1]
##  DelayedArray           0.12.2     2020-01-06 [2]
##  DEoptimR               1.0-8      2016-11-19 [1]
##  desc                   1.2.0      2018-05-01 [2]
##  destiny              * 3.0.1      2020-01-16 [1]
##  devtools               2.2.1      2019-09-24 [2]
##  digest                 0.6.25     2020-02-23 [1]
##  dplyr                * 1.0.2      2020-08-18 [1]
##  e1071                  1.7-3      2019-11-26 [1]
##  ellipsis               0.3.1      2020-05-15 [1]
##  evaluate               0.14       2019-05-28 [2]
##  fansi                  0.4.1      2020-01-08 [2]
##  farver                 2.0.3      2020-01-16 [1]
##  fitdistrplus           1.0-14     2019-01-23 [2]
##  forcats              * 0.5.0      2020-03-01 [1]
##  foreign                0.8-74     2019-12-26 [3]
##  formattable            0.2.0.1    2016-08-05 [1]
##  fs                     1.5.0      2020-07-31 [1]
##  future                 1.15.1     2019-11-25 [2]
##  future.apply           1.4.0      2020-01-07 [2]
##  gbRd                   0.4-11     2012-10-01 [2]
##  gdata                  2.18.0     2017-06-06 [2]
##  generics               0.0.2      2018-11-29 [2]
##  GenomeInfoDb           1.22.0     2019-10-29 [2]
##  GenomeInfoDbData       1.2.2      2020-01-14 [2]
##  GenomicRanges          1.38.0     2019-10-29 [2]
##  ggplot.multistats      1.0.0      2019-10-28 [1]
##  ggplot2              * 3.3.2      2020-06-19 [1]
##  ggrepel                0.8.1      2019-05-07 [2]
##  ggridges               0.5.2      2020-01-12 [2]
##  ggthemes               4.2.0      2019-05-13 [1]
##  globals                0.12.5     2019-12-07 [2]
##  glue                   1.3.2      2020-03-12 [1]
##  gplots                 3.0.1.2    2020-01-11 [2]
##  gridExtra              2.3        2017-09-09 [2]
##  gtable                 0.3.0      2019-03-25 [2]
##  gtools                 3.8.1      2018-06-26 [2]
##  haven                  2.3.1      2020-06-01 [1]
##  hexbin                 1.28.0     2019-11-11 [2]
##  hms                    0.5.3      2020-01-08 [1]
##  htmltools              0.4.0      2019-10-04 [2]
##  htmlwidgets            1.5.1      2019-10-08 [2]
##  httr                   1.4.2      2020-07-20 [1]
##  ica                    1.0-2      2018-05-24 [2]
##  igraph                 1.2.5      2020-03-19 [1]
##  IRanges                2.20.2     2020-01-13 [2]
##  irlba                  2.3.3      2019-02-05 [2]
##  jsonlite               1.7.1      2020-09-07 [1]
##  KernSmooth             2.23-16    2019-10-15 [3]
##  knitr                  1.26       2019-11-12 [2]
##  knn.covertree          1.0        2019-10-28 [1]
##  labeling               0.3        2014-08-23 [2]
##  laeken                 0.5.1      2020-02-05 [1]
##  lattice                0.20-38    2018-11-04 [3]
##  lazyeval               0.2.2      2019-03-15 [2]
##  leiden                 0.3.1      2019-07-23 [2]
##  lifecycle              0.2.0      2020-03-06 [1]
##  listenv                0.8.0      2019-12-05 [2]
##  lmtest                 0.9-37     2019-04-30 [2]
##  lsei                   1.2-0      2017-10-23 [2]
##  lubridate              1.7.9.2    2020-11-13 [1]
##  magrittr             * 2.0.1      2020-11-17 [1]
##  MASS                   7.3-51.5   2019-12-20 [3]
##  Matrix                 1.2-18     2019-11-27 [3]
##  matrixStats            0.56.0     2020-03-13 [1]
##  memoise                1.1.0      2017-04-21 [2]
##  metap                  1.2        2019-12-08 [2]
##  mnormt                 1.5-5      2016-10-15 [2]
##  modelr                 0.1.8      2020-05-19 [1]
##  multcomp               1.4-12     2020-01-10 [2]
##  multtest               2.42.0     2019-10-29 [2]
##  munsell                0.5.0      2018-06-12 [2]
##  mutoss                 0.1-12     2017-12-04 [2]
##  mvtnorm                1.0-12     2020-01-09 [2]
##  nlme                   3.1-143    2019-12-10 [3]
##  nnet                   7.3-12     2016-02-02 [3]
##  npsurv                 0.4-0      2017-10-14 [2]
##  numDeriv               2016.8-1.1 2019-06-06 [2]
##  openxlsx               4.1.5      2020-05-06 [1]
##  pbapply                1.4-2      2019-08-31 [2]
##  pcaMethods             1.78.0     2019-10-29 [2]
##  pillar                 1.4.6      2020-07-10 [1]
##  pkgbuild               1.0.6      2019-10-09 [2]
##  pkgconfig              2.0.3      2019-09-22 [1]
##  pkgload                1.0.2      2018-10-29 [2]
##  plotly                 4.9.1      2019-11-07 [2]
##  plotrix                3.7-7      2019-12-05 [2]
##  plyr                   1.8.5      2019-12-10 [2]
##  png                    0.1-7      2013-12-03 [2]
##  prettyunits            1.1.1      2020-01-24 [1]
##  processx               3.4.2      2020-02-09 [1]
##  progeny              * 1.11.3     2020-10-22 [1]
##  proxy                  0.4-24     2020-04-25 [1]
##  ps                     1.3.2      2020-02-13 [1]
##  purrr                * 0.3.4      2020-04-17 [1]
##  R.methodsS3            1.7.1      2016-02-16 [2]
##  R.oo                   1.23.0     2019-11-03 [2]
##  R.utils                2.9.2      2019-12-08 [2]
##  R6                     2.4.1      2019-11-12 [1]
##  ranger                 0.12.1     2020-01-10 [1]
##  RANN                   2.6.1      2019-01-08 [2]
##  rappdirs               0.3.1      2016-03-28 [2]
##  RColorBrewer           1.1-2      2014-12-07 [2]
##  Rcpp                   1.0.4      2020-03-17 [1]
##  RcppAnnoy              0.0.16     2020-03-08 [1]
##  RcppEigen              0.3.3.7.0  2019-11-16 [2]
##  RcppHNSW               0.2.0      2019-09-20 [2]
##  RcppParallel           4.4.4      2019-09-27 [2]
##  RCurl                  1.98-1.1   2020-01-19 [1]
##  Rdpack                 0.11-1     2019-12-14 [2]
##  readr                * 1.4.0      2020-10-05 [1]
##  readxl               * 1.3.1      2019-03-13 [2]
##  rematch                1.0.1      2016-04-21 [2]
##  remotes                2.1.0      2019-06-24 [2]
##  reprex                 0.3.0      2019-05-16 [2]
##  reshape2               1.4.3      2017-12-11 [2]
##  reticulate             1.14       2019-12-17 [2]
##  rio                    0.5.16     2018-11-26 [1]
##  rlang                  0.4.8      2020-10-08 [1]
##  rmarkdown              2.0        2019-12-12 [2]
##  robustbase             0.93-6     2020-03-23 [1]
##  ROCR                   1.0-7      2015-03-26 [2]
##  rprojroot              1.3-2      2018-01-03 [2]
##  RSpectra               0.16-0     2019-12-01 [2]
##  rstudioapi             0.11       2020-02-07 [1]
##  rsvd                   1.0.3      2020-02-17 [1]
##  Rtsne                  0.15       2018-11-10 [2]
##  rvest                  0.3.6      2020-07-25 [1]
##  S4Vectors              0.24.2     2020-01-13 [2]
##  sandwich               2.5-1      2019-04-06 [2]
##  scales                 1.1.0      2019-11-18 [2]
##  scatterplot3d          0.3-41     2018-03-14 [1]
##  sctransform            0.2.1      2019-12-17 [2]
##  SDMTools               1.1-221.2  2019-11-30 [2]
##  sessioninfo            1.1.1      2018-11-05 [2]
##  Seurat               * 3.1.2      2019-12-12 [2]
##  SingleCellExperiment   1.8.0      2019-10-29 [2]
##  smoother               1.1        2015-04-16 [1]
##  sn                     1.5-4      2019-05-14 [2]
##  sp                     1.4-2      2020-05-20 [1]
##  stringi                1.5.3      2020-09-09 [1]
##  stringr              * 1.4.0      2019-02-10 [1]
##  SummarizedExperiment   1.16.1     2019-12-19 [2]
##  survival               3.1-8      2019-12-03 [3]
##  testthat               2.3.2      2020-03-02 [1]
##  TFisher                0.2.0      2018-03-21 [2]
##  TH.data                1.0-10     2019-01-21 [2]
##  tibble               * 3.0.4      2020-10-12 [1]
##  tidyr                * 1.1.2      2020-08-27 [1]
##  tidyselect             1.1.0      2020-05-11 [1]
##  tidyverse            * 1.3.0      2019-11-21 [2]
##  tsne                   0.1-3      2016-07-15 [2]
##  TTR                    0.23-6     2019-12-15 [1]
##  usethis                1.5.1      2019-07-04 [2]
##  uwot                   0.1.5      2019-12-04 [2]
##  vcd                    1.4-7      2020-04-02 [1]
##  vctrs                  0.3.5      2020-11-17 [1]
##  VIM                    6.0.0      2020-05-08 [1]
##  viridis              * 0.5.1      2018-03-29 [2]
##  viridisLite          * 0.3.0      2018-02-01 [2]
##  withr                  2.3.0      2020-09-22 [1]
##  xfun                   0.12       2020-01-13 [2]
##  xml2                   1.3.2      2020-04-23 [1]
##  xts                    0.12-0     2020-01-19 [1]
##  XVector                0.26.0     2019-10-29 [2]
##  yaml                   2.2.1      2020-02-01 [1]
##  zip                    2.0.4      2019-09-01 [1]
##  zlibbioc               1.32.0     2019-10-29 [2]
##  zoo                    1.8-7      2020-01-10 [2]
##  source                                 
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  Bioconductor                           
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Github (clauswilke/colorblindr@1ac3d4d)
##  R-Forge (R 3.6.2)                      
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  Bioconductor                           
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Github (saezlab/progeny@94bea1c)       
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.3)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
## 
## [1] /home/uhlitzf/R/lib
## [2] /usr/local/lib/R/site-library
## [3] /usr/local/lib/R/library